% $Id: nlss2.tex 204 2009-04-20 01:59:36Z jsibert $
%
% Author: David Fournier
% Copyright (c) 2008 Regents of the University of California
%

%meagnification=1400
%\def\mysection#1{{\noindent\bf\medskip #1\medskip}}
\def\wha{\widehat\alpha}
\def\wht{\widehat\theta}
\def\Nipjp{N_{i+1,j+1}}
\def\Nipm{N_{i+1,m}}
\def\Nij{N_{ij}}
\def\nij{n_{ij}}
\def\Ei{E_i}
\def\Nio{N_{i1}}
\def\Nim{N_{im}}
\def\Nnpj{N_{n+1,j}}
\def\nnpj{n_{n+1,j}}
\def\Nimm{N_{i,m-1}}
\def\Zij{Z_{ij}}
\def\Fij{F_{ij}}
\def\Zim{Z_{im}}
\def\Zimm{Z_{i,m-1}}
\def\qip{q_{i+1}}
\def\qi{q_i}
\def\ui{u_i}
\def\sj{s_j}
\def\deltai{\delta_i}
\def\di{d_i}
\def\epsij{\epsilon_{ij}}
\def\kappaij{\kappa_{ij}}
\def\epst{\epsilon_t}
\def\at{\alpha_t}
\def\atp{\alpha_{t+1}}
\def\gami{\gam_i}
\def\lami{\lambda_i}
\def\etai{\eta_i}
\def\Cijobs{C^{obs}_{ij}}
\def\Cij{C_{ij}}
\def\myeq{}
\mysection{Introduction}

An enduring characteristic of computer science and technology
appears to be that hardware develops more quickly than software,
so that it is good practice to devlop todays software for tomorrow's
hardware. Currently it appears that one of the next dvelopments 
in computer hardware will be the widespread availability of mult-core processors so that soon
processors with eight, sixteen, or more processing cores will be commonplace.
The obvvious software devlepment to take advantage of this hardware devlopment 
is parallel processing, that is computer programs which do there computaitons in parallel. However the devlopment of parallel programming software 
which could permit the user to write custom parllel programs the has been 
disappointing. There has been more success in developing parallel code for
particular ``generic'' applications such as linear algebra 
where the class of situations to be dealt  with is well understood in advance.
The lesson is that to extend the applicability of parallel processing 
one should pick a class of models or programming problems which is 
well understood and implement the parallel aspects of the programming in such as
way that they they are as transparent as possible to the end user.
In this paper we describe such a programming paradigm for 
nonlinear (and linear) state space models. The programming is 
implemented using the random effects software module from
the AD Model Builder software package.

\def\B{\boldsymbol{0}}
\def\Bzero{\boldsymbol{0}}
\def\Byt{\boldsymbol{y_t}}
\def\Bytwo{\boldsymbol{y_2}}
\def\Bytwo{\boldsymbol{Y_2}}
\def\Byo{\boldsymbol{y_1}}
\def\Byone{\boldsymbol{y_1}}
\def\Bytwo{\boldsymbol{y_2}}
\def\ByT{\boldsymbol{y_T}}
\def\BY{\boldsymbol{Y}}
\def\Bdt{\boldsymbol{d_t}}
\def\Bct{\boldsymbol{c_t}}
\def\BZt{\boldsymbol{Z_t}}
\def\BQt{\boldsymbol{Q_t}}
\def\BTt{\boldsymbol{T_t}}
\def\BHt{\boldsymbol{H_t}}
\def\Balphat{\boldsymbol{\alpha_t}}
\def\BalphaT{\boldsymbol{\alpha_T}}
\def\BalphaTm{\boldsymbol{\alpha_{T-1}}}
\def\Balphatm{\boldsymbol{\alpha_{t-1}}}
\def\BalphaTmtwo{\boldsymbol{\alpha_{T-2}}}
\def\Balphatp{\boldsymbol{\alpha_{t+1}}}
\def\Balpha{\boldsymbol{\alpha}}
\def\Balphatwo{\boldsymbol{\alpha_2}}
\def\Balphathree{\boldsymbol{\alpha_3}}
\def\Balphao{\boldsymbol{\alpha_1}}
\def\Balphaone{\boldsymbol{\alpha_1}}
\def\Balphaz{\boldsymbol{\alpha_0}}
\def\Baz{\mathbf{a}_0}
\def\BPz{\mathbf{P}_0}
\def\Betat{\boldsymbol{\eta_t}}
\def\Balphato{\boldsymbol{\alpha_{t-1}}}
\def\Bepst{\boldsymbol{\epsilon_t}}
\def\Bepstprime{\boldsymbol{\epsilon_{t^\prime}}}
\def\Bepsone{\boldsymbol{\epsilon_1}}
\def\Bepstwo{\boldsymbol{\epsilon_2}}
\def\BepsT{\boldsymbol{\epsilon_T}}
\def\Beps{\boldsymbol{\epsilon}}

\mysection{The linear state space model}
To introduce the ideas we begin with the linear state space model where the
calculations are very simple. We have adopted the  notation used in Harvey.
For $t=1,\ldots,N$ let $\Byt$ be an  $N$ dimensional vector of observations.
These observations are related to an $m$ dimensional vector $\Balphat$
known as the state vector via the measurement equation
\begin{equation}
\Byt=\BZt(\theta)\Balphat+\Bdt(\theta)+\Bepst\label{nlss:xx01}
\end{equation}
where $\BZt$ is an $N\times m$ matrix, $\Bdt$ is an  $N$ dimensional vector
and $\Bepst$ is an $N$ dimensional multivariate normal random vector 
with mean $\Bzero$ and $N\times N$ covariance matrix $\BHt$.
The  $\BZt$  etc. depend on a vector of parameters $\theta$.


In general the states $\Balphat$ are not observable. They are assumed to be generated by a first-order Markov process.
\begin{equation}
\Balphat=\BTt(\theta)\Balphato+\Bct(\theta)+\Betat\label{nlss:xx02}
\end{equation}
where $\BTt$ in an $m \times m$ matrix $\Bct$ is an n dimensional vector
and $\Betat$ is an n dimensional multivariate normal random vector 
with mean $\Bzero$ and $m\times m$ covariance matrix $\BQt$.
The relationship \ref{nlss:xx02} is referredd to as the transition equation.
It can be used to solve for the
$\Balphat$ in terms of the $\Bepstprime$ for $t^\prime \le t$,
so that we write
$\Balphat(\Beps)$.

To complete the model specification assume that the intial state vector
$\Balphaz$ has mean $\Baz$ and covariance matrix $\BPz$.
From the linearity of equations
\ref{nlss:xx01} and \ref{nlss:xx02}
%it follows that $Byo$ and $\Balphao$ are jointly multivariate 
%noramlly distributed, that is the logarithm of their joint probability density 
%$G_t(y_1,\alpha_1)$ is quadratic in $y_1$ and $\alpha_1$ so that 
%the probaiblity density function for $\alpha_1$ given a value $y_1$
%for $\Byo$ can be found by maximizing $G_1(y_1,\alpha_1)$ with respect
%to $\alpha_1$, and the Hessian at the maximximg value determines the
%covariance matrix for $\Balphao$.  Of course for a quadratic function 
%this maximinzation can be carred out uising purely algebraic manipultions
%which perhaps obscure what is really going on.
%Having obtained the proability distribution for $\Balphao$
%given that $\Byo=y_1$ we can continue. We can calculate 
%the probability distribution for $\Balphat$ using equation
%\ref{nlss:xx02}. Now we can calculate the joint distribtuion for
%$\Bytwo$ and $\Balphatwo$ and given a value for
% a value for $\Bytwo$ calcuate the conditonal distribution of
%$\Balphatwo$ via maximization as before.
%It is clear that this can be carried on for all values of $t$ and
%this sequential estimation  process is known as the Kalman filter.
%The Kalman filter has several drawbacks. The first of these is that it
%is a sequential process which does not lend itself to parallel processing.
%Also it does not use all the information in the data since it only
%looks forward one time step. Aslo while it has good properlties for the
%normal, linear state space model when non normal distributions and nonlinear
%state space models are considered its properties are suboptimal. 

Let $F_t(\Byt,\alpha;\theta)$
be the probability density for $\Byt$ given $\Balphat$.
\begin{equation}
 F_t(y_t,\alpha;\theta)=|\BHt|^{-\frac{1}{2}}
  \exp\Big(-\frac{1}{2}
    \big(\Byt-\BZt(\theta)\Balphat(\epsilon)-\Bdt(\theta)\big)^\prime \BHt^{-1} 
  \big (\Byt-\BZt(\theta)\Balphat(\epsilon)-\Bdt(\theta)\big)\Big) 
\end{equation}
and $\Phi_t(\Bepst;\theta)$ be the proability density for $\Bepst$.
\begin{equation}
 \Phi_t(\Bepst;\theta)=|\BQt|^{-\frac{1}{2}}
  \exp\big(-\frac{1}{2} \Bepst^\prime \BQt^{-1} \Bepst \big )
\end{equation}
Let $\BY=(\Byone,\Bytwo,\ldots,\ByT)$ be all the observations.
The marginal distribution $L(\BY,\theta)$ for $\BY$ is found by integrating over the
$\Bepst$
\begin{equation}
L(Y,\theta)=\int \prod_t F_t(y_t,\alpha(\Beps);\theta)
 \Phi_t(\Bepst;\theta)\,d\Beps 
\end{equation}
It is more convenient to take logarithms to
get
\begin{equation}
L(\BY,\theta)=
  \int \exp\Big\{-\sum_t -\log\big(F_t(y_t,\alpha(\Beps);\theta)\big)
 -\log\big(\Phi_t(\Bepst;\theta)\big)\Big\}\,d\Beps \label{nlss:xx02a}
\end{equation}
The maximum likelihood estimates for 
$\wht$ are found by maximizing 
$L(\BY,\theta)$ with respect to~$\theta$.
\begin{equation}
\wht=\max_\theta \{ L(\BY,\theta)\}\label{nlss:xx03}
\end{equation}
The Kalman filter (Harvey) is a method for finding an approximate solution
to \ref{nlss:xx03}. While the Kalman filter has certain computational 
advantages it is inherently sequential and does not lend itself to 
parallelization. Also its properties degrade rapidily as nonlinear aspects 
are introduced into the model. In contrast the holistic method
provides the exact solution to \ref{nlss:xx03} for linear
models, is easily parallelizeable ,
and performs better for nonlinear models.
From the definition of the linear state space model it follows that
the expression
\begin{equation}
  \sum_t -\log\big(F_t(y_t,\alpha(\Beps);\theta)\big)
 -\log\big(\Phi_t(\Bepst;\theta)\big) \label{nlss:xx04}
\end{equation}
is a quadratic function of $\Beps$, so that the
integral can be calculated exactly by finding $\widehat\Beps$
which minimizes \ref{nlss:xx04} togehter with the Hessian 
with repsect to $\Beps$.
This procedure is know as the Laplace approximation to the integral
and for quadratic functions it is exact.
If we make the change of variables
\begin{equation}
 \Bepst=\Balphat -\BTt(\theta)\Balphato-\Bct(\theta)\label{nlss:xx05}
\end{equation}
the integral becomes
\begin{equation}
L(\BY,\theta)=
  \int \exp\Big\{-\sum_t -\log\big(F_t(y_t,\alpha_t;\theta)\big)
 -\log\big(\Phi_t(
   \Balphat -\BTt(\theta)\Balphato-\Bct(\theta);
     \theta)\big)\Big\}\,d\Balpha \label{nlss:xx06}
\end{equation}
In general a change of variables like \ref{nlss:xx05}
involves the matrix of the Jacobian of the transformation,
but due to the simple form for the linear state space model the
Jacobian matrix is equal to 1.
Note that \ref{nlss:xx05}
is the sum of $T$ terms where the $t$'th term involves only
$\Balphat$ and $\Balphato$. Each of these terms can be calculated
independently of the other ones. This is the basis of the
parallelizability of the method.


\mysection{The nonlinear state space model}
For the nonlinear state space model replace the transition equation
\ref{nlss:xx02}by the more general nonlinear relationship
\begin{equation}
 \myeq{\Balphat=g_t(\Balphato,\Bepst;\theta)}\label{nlss:xx1}
\end{equation}
where $\Bepst$
is an  random vector with probability density function
$\Phi_t(\Bepst;\theta)$. 
We do not require the $\Bepst$ be multivariate normal.
We also do not assume that the random vecotrs $\Balphat$
and $\Bepst$ have the same dimension for each value of $t$.

We replace the measurement equation \ref{nlss:xx01}
by the assumption that the 
$\Byt$ is a random vector with
proability density function $F_t(\Byt,\Balphat,\Balphatm;\theta)$.
It is convenicen for model formulaiton to allow $F_t$ to depend on
both $\Balphat$ and $\Balphatm$.
Otherwise we make no assumptions
about the particular form of $F_t$,
in particular no assumptions about normality.
Under these assumptions 
the marginal distribution $L(\BY,\theta)$ for $\BY$ is given by 
%\begin{equation}
\ref{nlss:xx02a} however 
\ref{nlss:xx04} is no longer quadratic  so that its minimization and the subsequent calculations are more difficult.

%However first 
%we will derive the equations for the Kalman filter for linear state space models and show how that connects via optimziation over $\alpha$ to our
%general approach.
%\mysection{The Kalman filter for linear state space models}
%\begin{equation}
% \myeq{
%  (y-Z\alpha-d)^\prime H^{-1}(y-Z\alpha-d)+(\alpha-e)^\prime S^{-1}(\alpha-e)
%}\label{nlss:xx100}
%\end{equation}
%Consider this as a quadratic function in $\alpha$.
%\begin{equation}
% \myeq{
%  h(\alpha)=A+ B\alpha +{1\over2}\alpha^{\prime}C\alpha 
%}\label{nlss:xx101}
%\end{equation}
%Now let $\wha=-C^{-1}B$, so that
%\begin{equation}
% \myeq{
%  h(\wha+\delta)=A+ B\wha +{1\over2}\wha^{\prime}C\wha 
%     +{1\over2}\delta^{\prime}C\delta
%}\label{nlss:xx102}
%\end{equation}
%which simplifies to
%\begin{equation}
% \myeq{
%  h(\wha+\delta)=A -{1\over2}\wha^{\prime}C\wha 
%     +{1\over2}\delta^{\prime}C\delta
%}\label{nlss:xx102}
%\end{equation}
%so the integral becomes
%\begin{equation}
% \myeq{
%   \int \exp( -A +{1\over2}\wha^{\prime}C\wha) 
%     \exp(-{1\over2}\delta^{\prime}C\delta)\,d\delta
%}\label{nlss:xx103}
%\end{equation}
%or
%\begin{equation}
% \myeq{
%   |C|^{-{1\over2}}\exp( -A +{1\over2}\wha^{\prime}C\wha) 
%}\label{nlss:xx103}
%\end{equation}
%since $\wha=-C^{-1}B$ this simplifies to 
%\begin{equation}
% \myeq{
%   |C|^{-{1\over2}}\exp( -A +{1\over2}B^{\prime}C^{-1}B) 
%}\label{nlss:xx103}
%\end{equation}
%solving for $A$ and $B$ we get
%\begin{equation}
% \myeq{
%   A=(y-d)^\prime H^{-1} (y-d) +e^\prime S^{-1} e
%}\label{nlss:xx103}
%\end{equation}
%and
%\begin{equation}
% \myeq{
% B=(y-d)^\prime H^{-1}Z+e^\prime S^{-1}
%}\label{nlss:xx104}
%\end{equation}
%and
%\begin{equation}
% \myeq{
% C=Z^\prime H^{-1}Z+S^{-1}
%}\label{nlss:xx105}
%\end{equation}
%and
%\begin{equation}
% \myeq{
% B^\prime C^{-1}B=\Big((y-d)^\prime H^{-1}Z+e^\prime S^{-1}\Big)\Big\{Z^\prime H^{-1}Z+S^{-1}\Big\}^{-1}
%  \Big(Z^\prime H^{-1}(y-d) +S^{-1}e\Big)
%}\label{nlss:xx106}
%\end{equation}
%
We replace the relationship
\ref{nlss:xx05} by the following assumption.
For any value of the parameters $\Balphato$
assume that it is possible to solve the  equations\ref{nlss:xx1}
for  $\Bepst$ in terms of $\Balphat$ ie
\begin{equation}
 \myeq{\Bepst=h_t(\Balphat,\Balphato)}\label{nlss:xx3}
\end{equation}
so that the map $\Balpha \rightarrow \Beps$  
is a diffeomorphism. The Jacobian matrix $J(\Balpha)$ of this map
is a banded upper triangular matrix of the form

\[
\left | \begin{array}{ccccc}
   D_1h_1(\Balphaone,\Balphaz) & 0  & 0 &  & 0 \\
   D_2h_2(\Balphatwo,\Balphaone) & D_1h_2(\Balphatwo,\Balphaone)  & 0  & & 0  \\ 
     0                    &  D_2h_3(\Balphathree,\Balphatwo)    &  D_1h_3(\Balphathree,\Balphatwo)  & & 0\\
                          &                           &                & \ddots \\
     0                    & 0                         &    0           &  D_1h_{T-1}(\BalphaTm,\BalphaTmtwo)           & 0 \\
     0                    & 0                         &    0           &  D_2h_T(\BalphaT,\BalphaTm)  & D_1h_T(\BalphaT,\BalphaTm) 
 
   \end{array}
  \right |
\]


% \[
% \left | \begin{array}{ccccc}
%    D_1h_1(\Balphatwo,\Balphaone) & D_2h_1(\alpha_2,\alpha_1)  & 0 &  & 0 \\
%      0                    & D_1h_2(\alpha_3,\alpha_2)  &  D_2h_2(\alpha_3,\alpha_2) & & 0  \\ 
%      0                    &  0                        &  D_1h_3(\alpha_4,\alpha_3)  & & 0\\
%                           &                           &                & \ddots \\
%      0                    & 0                         &    0           &           & D_2h_{n-1}(\alpha_n,\alpha_{n-1}) \\
%      0                    & 0                         &    0           &           & D_1h_n(\alpha_{n+1},\alpha_n) 
%  
%    \end{array}
%   \right |
% \]
% 
where $D_1h_t(\Balphat,\Balphato)$ is the (square) matrix of partial derivatvies of $h_t$ with respect to  $\Balphat$ and 
 $D_2h_t(\Balphat,\Balphato)$ is the matrix of partial derivatives 
of $h_t$ with respect to  $\Balphato$. 
We can use the $h_t$  to
repararameterize the integral  in terms of the $\alpha_t$. 
This reparameterization involves the determinant of the Jacobian matrix, $J$.
 
\begin{equation}
 \myeq{
    \int \prod_t F_t(\Byt,\Balpha)
    |J(\Balpha)|^{-1}\prod_t \Phi_t(h_t\big(\Balphat,\Balphatm)\big) d\Balpha}\label{nlss:xx5}
\end{equation}
\mysection{Approximate integration via the Laplace approximation}
  The laplace approximation for the integral \ref{nlss:xx5}
  involves approximating the integrand by the second order Taylor expansion 
  in $\alpha$ evaluated at the maximzing value $\widehat \alpha$.
    It is more convenient to maximize the log of the integrand or
   equivalently to minimize minus the log of the integrand.
    Let $f_t(\Byt,\Balphat)=\ln\big(F_t(\Byt,\Balphat)\big)$
    and $\phi_t(h_t\big(\Balphat,\Balphatm)=
            \ln\big(\Phi_t(h_t\big(\Balphat,\Balphatm)\big)$
We can write the integral as
\begin{equation}
 \myeq{
    \int \exp\Big( -\left\{-\sum_t f_t(y_t,\alpha_t)
    +\ln\big(|J(\alpha)|\big)
     -\sum_t \phi_t(h_t\big(\alpha_t,\alpha_{t-1})\big) \right\}\Big)
    \,d\alpha}
    \label{nlss:xx6}
\end{equation}
The computationally difficult term in \ref{nlss:xx6} is 
    $\ln\big(|J(\alpha)|\big)$.
Since $J(\alpha)$ is upper diagonal this is given by 
\begin{equation}
 \myeq{
    \ln\big(|J(\alpha)|\big)=\sum_t \ln(|D_1h_t(\alpha_t,\alpha_{t-1}|)}
    \label{nlss:xx7}
\end{equation}
which only involes the diagonal terms so that is it is not necessary to 
calculate the terms $|D_2h_t(\alpha_t,\alpha_{t-1})|$.
Gathering together terms with the sam $t$ subscript we obtain
\begin{equation}
 \myeq{
    \int \exp\Big( -\left\{-\sum_t f_t(y_t,\alpha_t)
     +\ln(|D_1h_t(\alpha_t,\alpha_{t-1})|
     -\phi_t(h_t\big(\alpha_t,\alpha_{t-1})\big) \right\}\big)d\alpha}
    \label{nlss:xx8}
\end{equation}
Let $\psi_t$ be defined by
\begin{equation}
 \myeq{
   \psi_t(\alpha_t,\alpha_{t-1};\theta)=-f_t(y_t,\alpha_t;\theta)
     +\ln\big (|D_1h_t(\alpha_t,\alpha_{t-1};\theta)|\big)
     -\phi_t\big (h_t(\alpha_t,\alpha_{t-1};\theta)\big)}
    \label{nlss:xx9}
\end{equation}
and rewrite the integral in terms of the $\psi_t$
\begin{equation}
 \myeq{
    \int 
      \exp\big( -\bigg\{-\sum_t \psi_t(\alpha_t,\alpha_{t-1};\theta)\bigg\}\big)
    \,d\alpha}
    \label{nlss:xx10}
\end{equation}
Now setting 
\begin{equation}
 \myeq{
   \Psi(\alpha;\theta)=\sum_t \psi_t(\alpha_t,\alpha_{t-1};\theta)}
    \label{nlss:xx11}
\end{equation}
so that the integral becomes
\begin{equation}
 \myeq{
    \int \exp\big( -\Big\{-\Psi(\alpha;\theta)\Big\}\big)\,d\alpha}
    \label{nlss:xx12}
\end{equation}
Let $\widehat\alpha$ be the value of $\alpha$ which minimizes
    $\Psi(\alpha;\theta)$ and approximate $\Psi(\alpha;\theta)$ 
 by its second order Taylor expansion at 
$\widehat\alpha$ we obtain
\begin{equation}
 \myeq{
    \int \exp\big( -\Big\{-\Psi(\widehat\alpha;\theta) 
     - D^2_\alpha\Psi(\wha;\theta)\Big\}\big)\,d\alpha}
    \label{nlss:xx13}
\end{equation}
which is equal to 
\begin{equation}
\exp(\Psi(\wha;\theta)|-D^2_\alpha\Psi(\wha;\theta)|^{1/2}
\end{equation}
where $|-D^2\Psi(\wha;\theta)|$ is the 
determinant of the Hessian of $-\Psi(\wha;\theta)$.

$\Psi$ is separable function of $\Balpha$ that is it is
the sum of $n$ terms each of which involves at most two of the $\alpha_t$.
This implies that the Hessian 
matrix is   block banded.  The block banded structure can be 
exploited in two ways. Calculating the determinant of a block diagonal matrix
is very efficient. Also solving the system of equations 
$u=H^{-1}v$ where $H$ is a positive definite symmetric block diagonal matrix
can be done very efficiently. This latter property can be exploited to
help solve for $\wha$ via the newton-raphon method and its variants.
 
These techniques are incorporated into the AD Model Builder
random effects software module. It is only necessary for the user to
formulate the problem in the correct fashion using the 
{\tt SEPARABLE\_FUNCTION} option. 

\beginexampledf
DATA_SECTION
  init_int n
  init_matrix Y(0,n,1,3)
PARAMETER_SECTION
  init_vector alpha0(1,3)
  init_bounded_number log_sigma_e(-2.0,5.0)
  init_bounded_number log_sigma_y(-2.0,5.0)
  random_effects_matrix alpha(1,n,1,3)
  matrix eps(1,n,1,3)
  objective_function_value f
PROCEDURE_SECTION
  f0(alpha0,log_sigma_y);
  int i;
  f1(alpha(1),alpha0,log_sigma_y,log_sigma_e,1);
  for (i=2;i<=n;i++)
  {
    f1(alpha(i),alpha(i-1),log_sigma_y,log_sigma_e,i);
  }

SEPARABLE_FUNCTION void f0(const dvar_vector& alpha0,const prevariable& log_sigma_y)
  dvariable s_y=exp(log_sigma_y);
  f+=3.0*log_sigma_y+0.5*norm2((Y(0)-alpha0)/s_y);
  
SEPARABLE_FUNCTION void f1(const dvar_vector& alphai,const dvar_vector& alphai1,const prevariable& log_sigma_y,const prevariable log_sigma_e,int i)
  dvariable s_e=exp(log_sigma_e);
  dvariable s_y=exp(log_sigma_y);
  dvar_vector epsi=alphai-alphai1;
  f+=3.0*log_sigma_e+0.5*norm2(epsi/s_e);
  f+=3.0*log_sigma_y+0.5*norm2((Y(i)-alphai)/s_y);

\endexampledf

% \mysection{Advantages of the holistic approach}
% 
% For the linear model with normally distributed random variables
% these are the full information maximum likelihood estimates. 
% 
% For nonlinear models and/or with non-normal random variables
%  the holistic approach produces better estimators than the Kalman filter.
% 
% The holistic approach is easily parallelizable so it 
% is a  good candidate for a computational paradigm which can exploit
% future multi-core computers.
% 
% Consider expression \ref{nlss:xx10}. 
% To calculate the Laplace approximaiton it 
% is necessary to maximize this expression with respect to the parameters 
% $\alpha_t$. This is the inner optimization. It has a separable
% strcture that is it is a sum of  summands $\psi_t$
% where each summand depends only on
% $\alpha_t$, $\alpha_{t-1}$, and
% $\theta$. As functions of these parameters the calcuations for each $\psi_t$
% are independent of each other and so can be trivially parallelized,
% and depend on only a relatively small number of parameters.
% After (and optionally during) the inner optimziation it is 
% also necessary to calculate the banded Hessian matrix, and to
% caluclate its determinant. These calculations are also 
% parallelizable and constitute all the calculations necessary tocarry out 
% one calculation  of $L(\theta)$. Finally we need to calculate the
% derivatives of $L(\theta)$ with respect to the parameters $\theta$.
% These derivatives can also be calculated within the 
% parallelizing structure used for the inner optimziation,
% and share the reduction in computational complexity due to the separable
% structure in the model.
% 

\mysection{A simple age-structured model}
This is an example of a simple model used in fisheries
stock assessment. There are $n$ years of catch data and the 
population is divided up into $m$ age classes where the last age 
class consists of all individual of age $m+a_0$ and older, 
where $a_0$ is the age of the
individuals in the first age class.
First we formulate the model in a manner which follows the 
stanbdard state space approach. 
Let $\Nij$ for $1\le i\le n+1$ and $1\le j\le m$ be the number of age class 
N indiviudual at the beginning of year $i$. 
For $1\le i\le n$ and $1\le j\le m$ the instantaneous total mortality rate
$\Zij$ is related to the $\Nij$ by the relationships
\begin{equation}
 \myeq{
   \Nipjp=\Nij\exp(-\Zij) \qquad \hbox{for}\qquad 1\le i\le n,\quad 1\le j<m-1}
    \label{nlss:xx200}
\end{equation}
and since the last age class is all the older fish
\begin{equation}
 \myeq{
   \Nipm= \Nimm\exp(-\Zimm)+ \Nim\exp(-\Zim)\qquad \hbox{for}\qquad 1\le i\le n}
    \label{nlss:xx201}
\end{equation}
Now the conventional way to parameterize the model is to start with
the number of fish at the beginning of year $1$ 
and the number of age class 1 fish at the beginning of each
$\Nio$ (recruitment) 
and to use the 
relationships \ref{nlss:xx200} and \ref{nlss:xx201} to compute
the other $\Nij$. 
It is of course necessary to supply values for the $Zij$, which we will deal 
with below, but the point is that one
starts with the fish population in a cetain "state" at the beginning of year
$1$ and calculates the later states as a function of this initial state and
its transitions from one year to the next. Now assume that
the total mortltiy can be decomposed into   component due to fishing, $\Fij$
and a component due to all other factors, $M$. 
where $\Fij$ is referred to as 
the instantaneous fishing mortality rate
and $M$ is referred to as the
the instantaneous natural mortality rate.
It is assumed that
\begin{equation}
 \myeq{
   \Zij=\Fij+M 
   }
    \label{nlss:xx202}
\end{equation}
It is further assumed that
\begin{equation}
 \myeq{
   \Fij = \qi \sj \di \Ei \exp(\kappaij) 
   }
    \label{nlss:xx203}
\end{equation}
where the $\Ei$ are exogenous parameters referred to as the fishing effort.
\begin{equation}
 \myeq{
   \qip=\qi\exp(\lami) \qquad \hbox{for} \quad 1<i<=n
  }
    \label{nlss:xx204}
\end{equation}
In addition assume the the recruitments $\Nio$ satify the relationship 
\begin{equation}
 \myeq{
   \Nio=R\exp(\deltai) \qquad \hbox{for} \quad 1<i<=n
  }
\end{equation}

\begin{equation}
 \myeq{
   \di=\exp(\etai) \qquad \hbox{for} \quad 1<i<=n
  }
    \label{nlss:xx205}
\end{equation}
The parameters $\etai$, $\kappaij$, $\deltai$, 
and $\lami$ are the random variables, 
that is, taken together they represent the
random variable $\Beps$ of the nonlinear state space model
which we wish to integrate over to
obtain the marginal distribution.

Let $\Cij$ be the catch of age class $j$ fish in year $i$.
which are given by the relationship.
\begin{equation}
 \myeq{
   \Cij={\Fij\over \Zij}\big(1-\exp(-\Zij)\big)\Nij
     \qquad \hbox{for} \quad 1\le i\le n \quad 1\le j\le m
  }
  \label{nlss:xx206}
\end{equation}
We also need some observational data to fit the model. We shall assume that
we have estimates of the catch at age $\Cijobs$.

\mysection{Holistic parameterization for the age-structured model}
The total number of random variables is
$nr$ $\kappaij$ parameters, $n$ $\etai$ parameters,
$n$ $\lami$ parameters, and $n-1$ $\deltai$ parameters,
so there are $nr+3n-1$ random variables in total.

We shall reparameterize the random variables by 
$N_{11}$, 
$\Nij$ for $2\le i\le n+1$, $1\le j\le m$, 
$\Nnpj$ for  $2\le j\le m$, 
$d_i$ for $1\le i\le n$,
$\qi$ for $2\le i\le n$, and $\ui$ for $1\le i\le n$.
where the $\ui$ are needed to split up the $m$'th age class.

\setcounter{MaxMatrixCols}{20}
\begin{equation}
\begin{matrix}
  & u_{11} & u_{12} & \ldots   & u_{1m} & u_{21} & u_{22} & \ldots & u_{2m} & u_{2,m+1} & u_{2,m+2} & u_{31} & u_{32} & u_{33} & \ldots & u_{3m}  \\
\noalign{\smallskip}
\delta_1 & 1 & 0  & \ldots & 0 & 0 \\
\chi_2 & 0 & 1 &  \ldots & 0 & 0 \\
\chi_3 & 0 & 0 & \ldots & 0 & 0 \\
\vdots \\
\chi_m & 0 & 0 & \ldots & 1  & 0 \\
\kappa_{11} &1/F_{11} & 0 & \ldots & 0 & -1/F_{11} & 0 \\
\kappa_{12} & 0 &1/F_{12} & \ldots & 0  & 0 &-1/F_{12} \\
\vdots \\
\kappa_{1m} & 0 & 0 & \ldots & 1/F_{1m}  & 0 & 0 &\ldots & -1/F_{1m}  \\
d_1 & 0 & 0 & \ldots & 0\\
u_1 & 0 & 0 & \ldots & 0 \\
\delta_1 & 0 & 0  & \ldots & 0 & 1\\

\end{matrix}
\end{equation}

